1 Loading the data

  1. Socio-demographic dta came from BRFSS 2000-2019.
  2. State-level CVD mortality rates were obtained from CDC WONDER.
  3. Data of primary care clinicians and cardiologists per 100,000 residents was  calculated from the Health Resources and Services Administration Area Health Resource.
  4. States’ political orientation.

1.1 BRFSS data

bfrss_overall <- read_rds(here('data', 'raw_data', 'bfrss_overall.RData')) %>% 
  mutate(state_id=as.numeric(state_id),
         year=as.numeric(year))

bfrss_race <- read_rds(here('data', 'raw_data', 'bfrss_race.RData')) %>% 
  mutate(state_id=as.numeric(state_id),
         year=as.numeric(year))

bfrss_gender <- read_rds(here('data', 'raw_data', 'bfrss_gender.RData'))%>% 
  mutate(state_id=as.numeric(state_id),
         year=as.numeric(year))

1.2 CDC WONDER data

cvd_overall <- read_tsv(here('data', 'raw_data', "state.txt")) %>% 
  select(c(State, "State Code",  "Year", "Age Adjusted Rate")) %>% 
  rename(state_id='State Code',
         year=Year,
         death='Age Adjusted Rate') %>% 
  mutate(state_id=as.numeric(state_id)) %>% 
  dplyr::filter(state_id !=11)#DC

cvd_black <- read_tsv(here('data', 'raw_data', "black.txt")) %>% 
  select(c(State, "State Code",  "Year", "Age Adjusted Rate")) %>% 
  rename(state_id='State Code',
         year=Year,
         death='Age Adjusted Rate')%>% 
  mutate(state_id=as.numeric(state_id),
         death=as.numeric(death))%>% 
  dplyr::filter(state_id !=11)#DC

cvd_white <- read_tsv(here('data', 'raw_data', "white.txt")) %>% 
  select(c(State, "State Code",  "Year", "Age Adjusted Rate")) %>% 
  rename(state_id='State Code',
         year=Year,
         death='Age Adjusted Rate')%>% 
  mutate(state_id=as.numeric(state_id),
         death=as.numeric(death))%>% 
  dplyr::filter(state_id !=11)#DC

cvd_hispanic <- read_tsv(here('data', 'raw_data', "hispanic.txt")) %>% 
  select(c(State, "State Code",  "Year", "Age Adjusted Rate")) %>% 
  rename(state_id='State Code',
         year=Year,
         death='Age Adjusted Rate')%>% 
  mutate(state_id=as.numeric(state_id),
         death=as.numeric(death))%>% 
  dplyr::filter(state_id !=11)#DC

cvd_men <- read_tsv(here('data', 'raw_data', "men.txt")) %>% 
  select(c(State, "State Code",  "Year", "Age Adjusted Rate")) %>% 
  rename(state_id='State Code',
         year=Year,
         death='Age Adjusted Rate')%>% 
  mutate(state_id=as.numeric(state_id),
         death=as.numeric(death))%>% 
  dplyr::filter(state_id !=11)#DC

cvd_women <- read_tsv(here('data', 'raw_data', "women.txt")) %>% 
  select(c(State, "State Code",  "Year", "Age Adjusted Rate")) %>% 
  rename(state_id='State Code',
         year=Year,
         death='Age Adjusted Rate')%>% 
  mutate(state_id=as.numeric(state_id),
         death=as.numeric(death))%>% 
  dplyr::filter(state_id !=11)#DC

1.3 Clinician data

clinician <- read_csv(here('data', 'raw_data', 'clinician.csv'))

1.4 States’ political orientation

politics <- read_excel(here('data', 'raw_data', 'political party.xlsx')) 
#0 republican, 1 democratic, 2 split

2 Cleaning and merging the data

2.1 Overall

overall <- list(cvd_overall, clinician, bfrss_overall) %>% 
  reduce(left_join, by = c('state_id', 'year'))%>% 
  filter(year>1999) %>% 
  left_join(politics, by=c('State', 'year')) %>% 
  mutate(#create exposure status
    post = case_when(State == "Alaska" & year < 2015 ~ 0,
                     State == "Idaho" & year < 2020 ~ 0,
                     State == "Indiana" & year<2015 ~ 0,
                     State == "Louisiana" & year<2016 ~ 0,
                     State == "Maine" & year<2019 ~ 0,
                     State == "Missouri" & year < 2020 ~ 0,
                     State == "Montana" & year<2016 ~ 0,
                     State == "Nebraska" & year<2020 ~ 0,
                     State == "Pennsylvania" & year<2015 ~ 0,
                     State == "Utah" & year < 2020 ~ 0,
                     State == "Virginia" & year < 2019 ~ 0,
                     !(State %in% c('Alaska', 'Idaho', "Indiana", 'Lousiana', "Maine", 'Missouri', "Montana", 
                                    "Nebraska", "Pennsylvania", "Utah", "Virginia")) & year<2014 ~ 0,
                     TRUE ~ 1), 
    treated = ifelse(State %in% c('Alabama', "Florida", "Georgia", "Mississippi", 
                                  "North Carolina", 'Oklahoma', "South Carolina", "South Dakota", 
                                  "Tennessee", "Texas", "Wisconsin", "Wyoming"), 0, 1),
    treatedpost = treated*post)
write_rds(overall, here("data", "output_data", "overall.rds"))

2.2 Black

black <- bfrss_race %>% 
  dplyr::filter(race=='Black') %>% 
  list(., clinician, cvd_black) %>% 
  reduce(right_join, by = c('state_id', 'year'))%>% 
  left_join(politics, by=c('State', 'year')) %>% 
  filter(year>1999) %>% 
  mutate(#create exposure status
    post = case_when(State == "Alaska" & year < 2015 ~ 0,
                     State == "Idaho" & year < 2020 ~ 0,
                     State == "Indiana" & year<2015 ~ 0,
                     State == "Louisiana" & year<2016 ~ 0,
                     State == "Maine" & year<2019 ~ 0,
                     State == "Missouri" & year < 2020 ~ 0,
                     State == "Montana" & year<2016 ~ 0,
                     State == "Nebraska" & year<2020 ~ 0,
                     State == "Pennsylvania" & year<2015 ~ 0,
                     State == "Utah" & year < 2020 ~ 0,
                     State == "Virginia" & year < 2019 ~ 0,
                     !(State %in% c('Alaska', 'Idaho', "Indiana", 'Lousiana', "Maine", 'Missouri', "Montana", 
                                    "Nebraska", "Pennsylvania", "Utah", "Virginia")) & year<2014 ~ 0,
                     TRUE ~ 1), 
    treated = ifelse(State %in% c('Alabama', "Florida", "Georgia", "Mississippi", 
                                  "North Carolina", 'Oklahoma', "South Carolina", "South Dakota", 
                                  "Tennessee", "Texas", "Wisconsin", "Wyoming"), 0, 1),
    treatedpost = treated*post) 
write_rds(black, here("data", "output_data", "black.rds"))

2.3 Hispanic

hispanic <- bfrss_race %>% 
  dplyr::filter(race=='Hispanic') %>% 
  list(., clinician, cvd_hispanic) %>% 
  reduce(right_join, by = c('state_id', 'year'))%>% 
  left_join(politics, by=c('State', 'year')) %>% 
  filter(year>1999) %>% 
  mutate(#create exposure status
    post = case_when(State == "Alaska" & year < 2015 ~ 0,
                     State == "Idaho" & year < 2020 ~ 0,
                     State == "Indiana" & year<2015 ~ 0,
                     State == "Louisiana" & year<2016 ~ 0,
                     State == "Maine" & year<2019 ~ 0,
                     State == "Missouri" & year < 2020 ~ 0,
                     State == "Montana" & year<2016 ~ 0,
                     State == "Nebraska" & year<2020 ~ 0,
                     State == "Pennsylvania" & year<2015 ~ 0,
                     State == "Utah" & year < 2020 ~ 0,
                     State == "Virginia" & year < 2019 ~ 0,
                     !(State %in% c('Alaska', 'Idaho', "Indiana", 'Lousiana', "Maine", 'Missouri', "Montana", 
                                    "Nebraska", "Pennsylvania", "Utah", "Virginia")) & year<2014 ~ 0,
                     TRUE ~ 1), 
    treated = ifelse(State %in% c('Alabama', "Florida", "Georgia", "Mississippi", 
                                  "North Carolina", 'Oklahoma', "South Carolina", "South Dakota", 
                                  "Tennessee", "Texas", "Wisconsin", "Wyoming"), 0, 1),
    treatedpost = treated*post) 

write_rds(hispanic, here("data", "output_data", "hispanic.rds"))

#Impute for Missouri, North Carolina, Utah, Wisconsin, Kansas, Maryland, Oregon, Illinois
hispanic_fill <- data.frame(State = c('Missouri','North Carolina', 'Utah','Wisconsin', 'Kansas', 'Maryland', 'Oregon',
                                      'Illinois'),
                            year=rep(2000:2019, each = 8)) %>% 
  merge(hispanic, by=c('State', 'year'), all.x=TRUE) %>% 
  fill(everything(), .direction='downup') 

hispanic_filled <- hispanic %>% filter(!State %in% c('Missouri','North Carolina', 'Utah','Wisconsin', 'Kansas', 
                                                    'Maryland', 'Oregon', 'Illinois')) %>% 
  bind_rows(hispanic_fill)
write_rds(hispanic_filled, here("data", "output_data", "hispanic_filled.rds"))

2.4 White

white <- bfrss_race %>% 
  dplyr::filter(race=='White') %>% 
  list(., clinician, cvd_white) %>% 
  reduce(right_join, by = c('state_id', 'year'))%>% 
  left_join(politics, by=c('State', 'year')) %>% 
  filter(year>1999) %>% 
  mutate(#create exposure status
    post = case_when(State == "Alaska" & year < 2015 ~ 0,
                     State == "Idaho" & year < 2020 ~ 0,
                     State == "Indiana" & year<2015 ~ 0,
                     State == "Louisiana" & year<2016 ~ 0,
                     State == "Maine" & year<2019 ~ 0,
                     State == "Missouri" & year < 2020 ~ 0,
                     State == "Montana" & year<2016 ~ 0,
                     State == "Nebraska" & year<2020 ~ 0,
                     State == "Pennsylvania" & year<2015 ~ 0,
                     State == "Utah" & year < 2020 ~ 0,
                     State == "Virginia" & year < 2019 ~ 0,
                     !(State %in% c('Alaska', 'Idaho', "Indiana", 'Lousiana', "Maine", 'Missouri', "Montana", 
                                    "Nebraska", "Pennsylvania", "Utah", "Virginia")) & year<2014 ~ 0,
                     TRUE ~ 1), 
    treated = ifelse(State %in% c('Alabama', "Florida", "Georgia", "Mississippi", 
                                  "North Carolina", 'Oklahoma', "South Carolina", "South Dakota", 
                                  "Tennessee", "Texas", "Wisconsin", "Wyoming"), 0, 1),
    treatedpost = treated*post) 

write_rds(white, here("data", "output_data", "white.rds"))

2.5 Men

men <- bfrss_gender %>% 
  dplyr::filter(male==1) %>% 
  list(., clinician, cvd_men) %>% 
  reduce(right_join, by = c('state_id', 'year'))%>% 
  filter(year>1999) %>% 
  left_join(politics, by=c('State', 'year')) %>% 
  mutate(#create exposure status
    post = case_when(State == "Alaska" & year < 2015 ~ 0,
                     State == "Idaho" & year < 2020 ~ 0,
                     State == "Indiana" & year<2015 ~ 0,
                     State == "Louisiana" & year<2016 ~ 0,
                     State == "Maine" & year<2019 ~ 0,
                     State == "Missouri" & year < 2020 ~ 0,
                     State == "Montana" & year<2016 ~ 0,
                     State == "Nebraska" & year<2020 ~ 0,
                     State == "Pennsylvania" & year<2015 ~ 0,
                     State == "Utah" & year < 2020 ~ 0,
                     State == "Virginia" & year < 2019 ~ 0,
                     !(State %in% c('Alaska', 'Idaho', "Indiana", 'Lousiana', "Maine", 'Missouri', "Montana", 
                                    "Nebraska", "Pennsylvania", "Utah", "Virginia")) & year<2014 ~ 0,
                     TRUE ~ 1), 
    treated = ifelse(State %in% c('Alabama', "Florida", "Georgia", "Mississippi", 
                                  "North Carolina", 'Oklahoma', "South Carolina", "South Dakota", 
                                  "Tennessee", "Texas", "Wisconsin", "Wyoming"), 0, 1),
    treatedpost = treated*post) 

write_rds(men, here("data", "output_data", "men.rds"))

2.6 Women

women <- bfrss_gender %>% 
  dplyr::filter(male==0) %>% 
  list(., clinician, cvd_women) %>% 
  reduce(right_join, by = c('state_id', 'year'))%>% 
  filter(year>1999) %>%
  left_join(politics, by=c('State', 'year')) %>% 
  mutate(#create exposure status
    post = case_when(State == "Alaska" & year < 2015 ~ 0,
                     State == "Idaho" & year < 2020 ~ 0,
                     State == "Indiana" & year<2015 ~ 0,
                     State == "Louisiana" & year<2016 ~ 0,
                     State == "Maine" & year<2019 ~ 0,
                     State == "Missouri" & year < 2020 ~ 0,
                     State == "Montana" & year<2016 ~ 0,
                     State == "Nebraska" & year<2020 ~ 0,
                     State == "Pennsylvania" & year<2015 ~ 0,
                     State == "Utah" & year < 2020 ~ 0,
                     State == "Virginia" & year < 2019 ~ 0,
                     !(State %in% c('Alaska', 'Idaho', "Indiana", 'Lousiana', "Maine", 'Missouri', "Montana", 
                                    "Nebraska", "Pennsylvania", "Utah", "Virginia")) & year<2014 ~ 0,
                     TRUE ~ 1), 
    treated = ifelse(State %in% c('Alabama', "Florida", "Georgia", "Mississippi", 
                                  "North Carolina", 'Oklahoma', "South Carolina", "South Dakota", 
                                  "Tennessee", "Texas", "Wisconsin", "Wyoming"), 0, 1),
    treatedpost = treated*post) 

write_rds(women, here("data", "output_data", "women.rds"))

3 Analyzing the data using the generalzied SCM

3.1 Overall

overall <- read_rds(here("data", "output_data", "overall.rds"))

est_overall <- gsynth(death ~ treatedpost + primary_rate + cardio_rate + 
                        population + 
              low_educ + employed_for_wages + party +  
              low_income + married + male + race_White,
            data = overall %>% filter(!state_id %in% c(2, 51)) %>% 
              filter(!(state_id==34 & year==2019)),  
            EM = F, 
            index = c("state_id","year"),
            inference = "parametric", se = TRUE, 
            #min.T0 = 5, 
            r = c(0, 5), 
            nboots = 200, CV = TRUE, force = "two-way", parallel = FALSE)

#ATT plot
p1 <- plot(est_overall, type = "counterfactual", raw = "none", 
           theme.bw = TRUE, main = "", #ylim = c(120, 200), 
           legendOff = TRUE, xlab = "", ylab = "") +
  theme(panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        #axis.text.x = element_text(angle = 90, vjust = 0.5),
        plot.margin = margin(10, 10, 50, 10))
#save plot
#ggsave(plot = p1, here("figures", 'overall att.png'))

overall_by_period <- est_overall[["est.att"]] %>% 
  as.data.frame() %>% 
  rownames_to_column(var= "time") %>% 
  mutate(time=as.numeric(time)) %>% 
  dplyr::filter(time>=0)
#write_csv(overall_by_period, here('tables', 'overall_by_period.csv'), row.names = FALSE) 


overall_average <- data.frame(est_overall$est.avg)
#write_csv(overall_average, here('tables', 'overall_average.csv'), row.names = FALSE)

3.2 Black

black <- read_rds(here("data", "output_data", "black.rds"))

est_black <- gsynth(death ~ treatedpost + primary_rate + cardio_rate + population + 
                      low_educ + employed_for_wages + party +  
                      low_income + married + male,
                      data = black %>% filter(!(State %in% c('Alaska', 'Hawaii','New Mexico', 
                                                            'Rhode Island', 'South Dakota', 'Utah', 'Virginia'))) %>% 
                                               filter(!(state_id==34 & year==2019)),  
                      EM = F, 
                      index = c("state_id","year"),
                      inference = "parametric", se = TRUE, 
                      #min.T0 = 5, 
                      r = c(0, 5), 
                      nboots = 200, CV = TRUE, force = "two-way", parallel = FALSE)

p2 <- plot(est_black, type = "counterfactual", raw = "none", theme.bw = TRUE, main = "",#ylim = c(220, 475),  
           legendOff = TRUE, xlab = "", ylab = "") +
  theme(panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        #axis.text.x = element_text(angle = 90, vjust = 0.5),
        plot.margin = margin(10, 10, 50, 10))
#save plot
#ggsave(plot = p2, here("figures", 'black att.png'))

black_by_period <- est_black[["est.att"]] %>% 
  as.data.frame() %>% 
  rownames_to_column(var= "time") %>% 
  mutate(time=as.numeric(time)) %>% 
  dplyr::filter(time>=0)

#write_csv(black_by_period, here('tables', 'black_by_period.csv'), row.names = FALSE) 

black_average <- data.frame(est_black$est.avg)

#write_csv(black_average, here('tables', 'black_average.csv'), row.names = FALSE)

3.3 Hispanic

hispanic <- read_rds(here("data", "output_data", "hispanic.rds"))
#gsynth doesn not work for Hispanic given the high missingness (so used the imputed data)

hispanic_filled <- read_rds(here("data", "output_data", "hispanic_filled.rds"))

est_hispanic <- gsynth(death ~ treatedpost + primary_rate + cardio_rate + population +
                       low_educ + employed_for_wages + party +
                       low_income + married + male,
                     data = hispanic_filled%>% filter(!State %in% c(
                       'Alabama', 'Arkansas', 'Idaho', 'Iowa', 'Kentucky',
                       'Louisiana', 'Minnesota', 'Nebraska', 'Rhode Island',
                       'South Carolina', 'Tennessee', 'Alaska', 'Delaware', 'Mississippi',
                       'New Hampshire', 'Virginia', 'Wyoming'
                       
                     )) %>% 
                       #%>% 
                       filter(!(state_id==34 & year==2019)),
                     EM = F,
                     index = c("state_id","year"),
                     inference = "parametric", se = TRUE,
                     #min.T0 = 5,
                     r = c(0, 5),
                     nboots = 200, CV = TRUE, force = "two-way", parallel = FALSE)

p3 <- plot(est_hispanic, type = "counterfactual", raw = "none", theme.bw = TRUE, main = "",#ylim = c(220, 475),
           legendOff = TRUE, xlab = "", ylab = "") +
   theme(panel.grid.major = element_blank(),
         panel.grid.minor = element_blank(),
         #axis.text.x = element_text(angle = 90, vjust = 0.5),
         plot.margin = margin(10, 10, 50, 10))
#save plot
#ggsave(plot = p3, here("figures", 'hispanic att.png'))

hispanic_by_period <- est_hispanic[["est.att"]] %>%
   as.data.frame() %>%
   rownames_to_column(var= "time") %>%
   mutate(time=as.numeric(time)) %>%
   dplyr::filter(time>=0)

#write_csv(hispanic_by_period, here('tables', 'hispanic_by_period.csv'), row.names = FALSE)

hispanic_average <- data.frame(est_hispanic$est.avg)

#write_csv(hispanic_average, here('tables', 'hispanic_average.csv'), row.names = FALSE)

3.4 White

white <- read_rds(here("data", "output_data", "white.rds"))

est_white <- gsynth(death ~ treatedpost + primary_rate + cardio_rate + population + 
                      low_educ + employed_for_wages + party +  
                      low_income + married + male,
                    data = white %>% filter(!state_id %in% c(2, 51))%>% 
                      filter(!(state_id==34 & year==2019)),  
                    EM = F, 
                    index = c("state_id","year"),
                    inference = "parametric", se = TRUE, 
                    #min.T0 = 5, 
                    r = c(0, 5), 
                    nboots = 200, CV = TRUE, force = "two-way", parallel = FALSE)

p4 <- plot(est_white, type = "counterfactual", raw = "none", theme.bw = TRUE, main = "",#ylim = c(220, 475),  
           legendOff = TRUE, xlab = "", ylab = "") +
  theme(panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        #axis.text.x = element_text(angle = 90, vjust = 0.5),
        plot.margin = margin(10, 10, 50, 10))
#save plot
#ggsave(plot = p4, here("figures", 'white att.png'))


white_by_period <- est_white[["est.att"]] %>% 
  as.data.frame() %>% 
  rownames_to_column(var= "time") %>% 
  mutate(time=as.numeric(time)) %>% 
  dplyr::filter(time>=0)

#write_csv(white_by_period, here('tables', 'white_by_period.csv'), row.names = FALSE) 

white_average <- data.frame(est_white$est.avg)

#write_csv(white_average, here('tables', 'white_average.csv'), row.names = FALSE)

3.5 Men

men <- read_rds(here("data", "output_data", "men.rds"))

est_men <- gsynth(death ~ treatedpost + primary_rate + cardio_rate + population + 
                      low_educ + employed_for_wages + party +  
                      low_income + married + race_White,
                    data = men %>% filter(!state_id %in% c(2, 51))%>% 
                    filter(!(state_id==34 & year==2019)),  
                    EM = F, 
                    index = c("state_id","year"),
                    inference = "parametric", se = TRUE, 
                    #min.T0 = 5, 
                    r = c(0, 5), 
                    nboots = 200, CV = TRUE, force = "two-way", parallel = FALSE)

p5 <- plot(est_men, type = "counterfactual", raw = "none", theme.bw = TRUE, main = "",#ylim = c(220, 475),  
           legendOff = TRUE, xlab = "", ylab = "") +
  theme(panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        #axis.text.x = element_text(angle = 90, vjust = 0.5),
        plot.margin = margin(10, 10, 50, 10))
#save plot
#ggsave(plot = p5, here("figures", 'men att.png'))

men_by_period <- est_men[["est.att"]] %>% 
  as.data.frame() %>% 
  rownames_to_column(var= "time") %>% 
  mutate(time=as.numeric(time)) %>% 
  dplyr::filter(time>=0)

#write_csv(men_by_period, here('tables', 'men_by_period.csv'), row.names = FALSE) 

men_average <- data.frame(est_men$est.avg)

#write_csv(men_average, here('tables', 'men_average.csv'), row.names = FALSE)

3.6 Women

women <- read_rds(here("data", "output_data", "women.rds"))

est_women <- gsynth(death ~ treatedpost + primary_rate + cardio_rate + population + 
                    low_educ + employed_for_wages + party +  
                    low_income + married + race_White,
                  data = women %>% filter(!state_id %in% c(2, 51))%>% 
                    filter(!(state_id==34 & year==2019)),  
                  EM = F, 
                  index = c("state_id","year"),
                  inference = "parametric", se = TRUE, 
                  #min.T0 = 5, 
                  r = c(0, 5), 
                  nboots = 200, CV = TRUE, force = "two-way", parallel = FALSE)

p6 <- plot(est_women, type = "counterfactual", raw = "none", theme.bw = TRUE, main = "",#ylim = c(220, 475),  
           legendOff = TRUE, xlab = "", ylab = "") +
  theme(panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        #axis.text.x = element_text(angle = 90, vjust = 0.5),
        plot.margin = margin(10, 10, 50, 10))
#save plot
#ggsave(plot = p6, here("figures", 'women att.png'))

women_by_period <- est_women[["est.att"]] %>% 
  as.data.frame() %>% 
  rownames_to_column(var= "time") %>% 
  mutate(time=as.numeric(time)) %>% 
  dplyr::filter(time>=0)
#write_csv(women_by_period, here('tables', 'women_by_period.csv'), row.names = FALSE) 

women_average <- data.frame(est_women$est.avg)
#write_csv(women_average, here('tables', 'women_average.csv'), row.names = FALSE)

5 Manuscript Figures and Tables

5.1 Table 1: Baseline characteristics

Baseline characheristics

overall <- read_rds(here("data", "output_data", "overall.rds"))
tab <- overall %>% 
filter(year==2014) %>% # keep the continuous var and the two categorical variables
  select(party, primary_rate, cardio_rate, 
         low_educ, employed_for_wages, party, death, 
         low_income, married, male, race_White, treated) %>% 
  tbl_summary(by = treated,
              missing = 'no',
              digits = list(all_continuous()~4,
                            all_categorical()~2), 
              statistic = list(all_continuous() ~ '{mean} ({sd})',
                               all_categorical() ~ '{n} / {N} ({p}%)')) %>% 
  bold_labels() %>% 
  #modify_header(stat_by = '**{level}**') %>% 
  modify_footnote(update = everything() ~ NA)

tab
Characteristic 0, N = 12 1, N = 38
party
    0 12.00 / 12.00 (100.00%) 12.00 / 38.00 (31.58%)
    1 0.00 / 12.00 (0.00%) 13.00 / 38.00 (34.21%)
    2 0.00 / 12.00 (0.00%) 13.00 / 38.00 (34.21%)
primary_rate 25.4078 (2.7933) 30.8041 (5.3187)
cardio_rate 5.3811 (1.4441) 6.6020 (2.3972)
low_educ 0.0895 (0.0352) 0.0609 (0.0265)
employed_for_wages 0.5970 (0.0738) 0.6524 (0.0621)
death 190.8333 (49.5080) 142.3105 (40.1301)
low_income 0.1152 (0.0321) 0.0904 (0.0279)
married 0.5988 (0.0390) 0.6190 (0.0449)
male 0.4099 (0.0309) 0.4282 (0.0274)
race_White 0.7291 (0.1180) 0.8043 (0.1291)
 
#gtsave(tab %>% as_gt(), here('tables', 'table1.png'), vwidth = 1250, vheight = 50)
#tab %>% as_flex_table() %>% flextable::save_as_docx(.,path=here("tables", "table1.docx"))

5.2 Table 2: Overall and subgroup effect

Overall and subgroup effect of the Medicaid expansion on CVD mortality

overall_average <- read_csv(here('tables', "overall_average.csv"))
overall_average$group='Overall'

black_average <- read_csv(here('tables', "black_average.csv"))
black_average$group='Black'

hispanic_average <- read_csv(here('tables', "hispanic_average.csv"))
hispanic_average$group='Hispanic'

white_average <- read_csv(here('tables', "white_average.csv"))
white_average$group='White'

men_average <- read_csv(here('tables', "men_average.csv"))
men_average$group='Men'

women_average <- read_csv(here('tables', "women_average.csv"))
women_average$group='Women'


average_data <- rbind(overall_average, black_average, hispanic_average, 
                 white_average, men_average, women_average)

average_data$group <- factor(average_data$group, 
                        levels=c('Overall', 'Black', 'Hispanic', 'White', 
                                 'Men','Women'))

average_table <- average_data %>% 
  transmute(Group = group,
            `Adjusted Mean Difference (95%CI)`= simulr::ci_fct(
              pe=Estimate,
              low=CI.lower,
              high=CI.upper,
              digit=2,
              multiplier = 1,
              scale = ""))

knitr::kable(average_table) %>%
  kableExtra::kable_styling(bootstrap_options = c("striped", "bordered"),
                full_width = F)
Group Adjusted Mean Difference (95%CI)
Overall -3.57 (-9.61, 2.46)
Black 5.24 (-10.45, 20.93)
Hispanic -6.89 (-23.51, 9.74)
White -3.88 (-8.15, 0.39)
Men -2.41 (-10.49, 5.67)
Women -3.87 (-7.61, -0.14) 
# average_table %>% flextable::flextable() %>% 
#   flextable::save_as_docx(.,path=here("tables", "table2.docx"))

5.3 Figure 1: Annual impact

Overall annual impact of Medicaid expansion on CVD deaths per 100,000 persons overall and by race and sex.

overall_by_period <- read_csv(here('tables', "overall_by_period.csv"))
black_by_period <- read_csv(here('tables', "black_by_period.csv"))
hispanic_by_period <- read_csv(here('tables', "hispanic_by_period.csv"))
white_by_period <- read_csv(here('tables', "white_by_period.csv"))
men_by_period <- read_csv(here('tables', "men_by_period.csv"))
women_by_period <- read_csv(here('tables', "women_by_period.csv"))

f_overall <- overall_by_period %>% 
  ggplot(aes(x = time, y = ATT)) +
  geom_point(color="#030303") + 
  geom_errorbar(aes( ymax = CI.lower, ymin = CI.upper) , width = .2, size = 0.7, 
                color="#030303")+
  geom_line(color="#030303") +
  ylim(-60, 60)+
  scale_x_continuous(breaks=seq(0,6,1))+
  geom_hline(yintercept = 0, lty=2,color="#030303") +
  ggtitle('Overall')+
  labs(x = "Number of years since Medicaid Expansion",
       y = "Change in CVD mortality, per 100,000 population") +
  theme_classic(base_size = 15) +
  theme(axis.line = element_line(colour = "grey50", size = 1),
        panel.background = element_blank(),
        axis.title = element_blank(),
        plot.title = element_text(hjust = 0.5))


f_black_hispanic_white <- bind_rows(black_by_period %>% mutate(group="black"),
                                    hispanic_by_period %>% mutate(group="hispanic"),
                                    white_by_period %>% mutate(group="white")) %>% 
  ggplot(aes(x = time, y = ATT, color = group)) +
  scale_color_manual(values=c("#8A2BE2", "#FF7F50", "#458B00"), 
                     labels=c('Black',   'Hispanic','White'))+
  geom_point(position = position_dodge(width = .5)) + 
  geom_errorbar( aes( ymax = CI.upper, ymin = CI.lower) , width = .2, 
                 position = position_dodge(width = .5), size = 0.7 )+
  geom_line(position = position_dodge(width = .5), size = 0.7) +
  ylim(-60, 60)+
  scale_x_continuous(breaks=seq(0,6,1))+
  #ggtitle('B')+
  geom_hline(yintercept = 0, lty=2, color="black") +
  theme(legend.position = "top",
        axis.title = element_blank(),
        axis.line = element_line(colour = "grey50", size = 1),
        panel.background = element_blank(),
        legend.background = element_rect(fill = "lemonchiffon", 
                                         colour = "grey50", 
                                         size = 1),
        legend.title = element_blank())

f_men_women <- bind_rows(men_by_period %>% mutate(group="men"),
                         women_by_period %>% mutate(group="women")) %>% 
  ggplot(aes(x = time, y = ATT, color = group)) +
  scale_color_manual(values=c("#1874CD", "#EE3B3B"), 
                     labels=c('Men','Women'))+
  geom_point(position = position_dodge(width = .5)) + 
  geom_errorbar( aes( ymax = CI.upper, ymin = CI.lower) , width = .2, 
                 position = position_dodge(width = .5), size = 0.7 )+
  geom_line(position = position_dodge(width = .5), size = 0.7) +
  ylim(-60, 60)+
  scale_x_continuous(breaks=seq(0,6,1))+
  #ggtitle('B')+
  geom_hline(yintercept = 0, lty=2, color="black") +
  theme(legend.position = "top",
        axis.title = element_blank(),
        axis.line = element_line(colour = "grey50", size = 1),
        panel.background = element_blank(),
        legend.background = element_rect(fill = "lemonchiffon", 
                                         colour = "grey50", 
                                         size = 1),
        legend.title = element_blank())

f_by_overall_race_sex <- ggarrange(f_overall, f_black_hispanic_white, f_men_women,
                                   ncol = 1, nrow=3) %>% 
  annotate_figure(left = textGrob("Change in CVD mortality, per 100,000 population", 
                                  rot = 90, vjust = 1, gp = gpar(cex = 1.3)), #cex = 0.8
                  bottom = textGrob("Number of years since Medicaid Expansion", gp = gpar(cex = 1.3)))

f_by_overall_race_sex

# ggsave(plot=f_by_overall_race_sex,
#        filename=here('figures', 'figure1_annual_impact_overall_race_sex.png'), 
#        width = 15, height = 25, units = "cm")

5.4 Figure 2. Pre-treatment fit

Pre-treatment fit showing the observed CVD deaths per 100,000 persons along with the counterfactual (synthetic control). Vertical line represents the beginning of expansion. The solid line is expansion state, dashed line is the synthetic control.

p1_6_fit <- ggarrange(
  p1+ rremove("ylab") + rremove("xlab")+ggtitle("Overall"), 
  p2+ rremove("ylab") + rremove("xlab")+ggtitle("Black"),
  p3+ rremove("ylab") + rremove("xlab")+ggtitle("Hispanic"),
  p4+ rremove("ylab") + rremove("xlab")+ggtitle("White"),
  p5+ rremove("ylab") + rremove("xlab")+ggtitle("Men"),
  p6+ rremove("ylab") + rremove("xlab")+ggtitle("Women"),
  font.label = list(size = 2, color = "black", face = "bold", family = NULL),
  ncol=2, nrow=3,common.legend = TRUE)


p1_6_fit

# ggsave(plot=p1_6_fit,
#        filename=here('figures', 'figure2_pre_treatment_fit.png'), 
#        width = 44, height = 36, units = "cm")

5.5 Figure 3: Analytical data structure

Analytical data structure of the Medicaid expansion states and control states for the overall study population.

q1 <- panelview(death ~treatedpost , data=overall, index=c("State","year"), 
                pre.post= TRUE,
                cex.axis.x=20,
                cex.axis.y=20,
                cex.lab=20)

q1

# ggsave(plot=q1, filename=here('figures', 'figure3_data_structure_overall.png'), 
#        width = 22, height = 18)

5.6 eTable 1: STROBE Statement—Checklist (See Appendix)

5.7 eTable 2: Numeric values for Figure 1

average_by_period_table <- 
  bind_rows(overall_by_period %>%  mutate(group="Overall"),
            black_by_period %>% mutate(group="Black"),
            hispanic_by_period %>% mutate(group="Hispanic"),
            white_by_period %>% mutate(group="White"),
            men_by_period %>% mutate(group="Men"),
            women_by_period %>% mutate(group="Women")) %>% 
  mutate(`Adjusted Mean Difference (95%CI)`= simulr::ci_fct(
    pe=ATT,
    low=CI.lower,
    high=CI.upper,
    digit=2,
    multiplier = 1,
    scale = "")) %>% 
  dplyr::select(group, time, `Adjusted Mean Difference (95%CI)`)
average_by_period_table %>% 
  knitr::kable() %>%
  kableExtra::kable_styling(bootstrap_options = c("striped", "bordered"),
                full_width = F)
group time Adjusted Mean Difference (95%CI)
Overall 0 -0.26 (-4.58, 4.07)
Overall 1 -6.05 (-12.35, 0.25)
Overall 2 -6.25 (-13.94, 1.45)
Overall 3 -1.30 (-10.74, 8.15)
Overall 4 -3.67 (-13.28, 5.94)
Overall 5 -3.26 (-13.84, 7.32)
Overall 6 -0.27 (-10.79, 10.25)
Black 0 7.20 (-7.06, 21.46)
Black 1 4.02 (-12.04, 20.08)
Black 2 -3.07 (-26.17, 20.03)
Black 3 4.62 (-22.23, 31.48)
Black 4 4.44 (-19.49, 28.38)
Black 5 8.65 (-14.57, 31.87)
Black 6 14.44 (-7.85, 36.74)
Hispanic 0 -6.78 (-18.63, 5.06)
Hispanic 1 -10.15 (-27.85, 7.56)
Hispanic 2 -8.36 (-30.07, 13.35)
Hispanic 3 -6.46 (-32.57, 19.65)
Hispanic 4 -1.98 (-27.05, 23.10)
Hispanic 5 -9.57 (-33.27, 14.12)
Hispanic 6 -4.43 (-29.24, 20.39)
White 0 -2.77 (-7.05, 1.50)
White 1 -5.86 (-10.99, -0.74)
White 2 -5.68 (-12.46, 1.10)
White 3 -2.43 (-9.96, 5.10)
White 4 -3.89 (-11.85, 4.07)
White 5 -3.24 (-9.86, 3.37)
White 6 -1.72 (-8.66, 5.22)
Men 0 0.74 (-5.71, 7.20)
Men 1 -5.01 (-15.06, 5.03)
Men 2 -6.60 (-18.09, 4.88)
Men 3 2.89 (-9.19, 14.97)
Men 4 -4.03 (-17.49, 9.43)
Men 5 -0.67 (-15.42, 14.08)
Men 6 -0.52 (-16.10, 15.05)
Women 0 -2.15 (-6.64, 2.33)
Women 1 -6.50 (-11.12, -1.87)
Women 2 -4.63 (-10.55, 1.29)
Women 3 -3.97 (-9.58, 1.63)
Women 4 -2.32 (-8.45, 3.81)
Women 5 -5.51 (-11.54, 0.52)
Women 6 0.37 (-5.61, 6.34) 
average_by_period_table %>% 
  pivot_wider(names_from = time, 
              values_from = `Adjusted Mean Difference (95%CI)`) %>% 
  flextable::flextable() %>% 
  flextable::save_as_docx(.,path=here("tables", "etable2.docx"))

5.8 eFigure 1: Data structure (black, hispanic, white)

Analytical data structure of Medicaid expansion states and control states for Blacks, Hispanics and Whites.

q2 <- panelview(death ~treatedpost , data=black, index=c("State","year"), 
                pre.post= TRUE)

#ggsave(here('figures', 'PV_black.png'), width = 15, height = 10)

q3 <- panelview(death ~treatedpost , data=hispanic, index=c("State","year"),
                pre.post= TRUE)

#ggsave(here('figures', 'PV_hispanic.png'), width = 15, height = 10)

q4 <- panelview(death ~treatedpost , data=white, index=c("State","year"), 
                pre.post= TRUE)

#ggsave(here('figures', 'PV_white.png'), width = 15, height = 10)

ggarrange(q2+ rremove("ylab") + rremove("xlab")+ggtitle("Black"), 
          q3+ rremove("ylab") + rremove("xlab")+ggtitle("Hispanic"),
          q4+ rremove("ylab") + rremove("xlab")+ggtitle("White"),
          ncol=1, nrow=3,common.legend = TRUE,
          cex.axis.x=20,
          cex.axis.y=20,
          cex.lab=20)
## $`1`

## 
## $`2`

## 
## attr(,"class")
## [1] "list"      "ggarrange"
# ggsave(here('figures', 'efigure1_data_structure_black_hispanic_white.png'), 
#        width = 36, height = 44, units = "cm")

5.9 eFigure 2: Data structure (men and women)

Analytical data structure of Medicaid expansion states and control states for men and women.

q5 <- panelview(death ~treatedpost , data=men, index=c("State","year"), pre.post= TRUE)

#ggsave(here('figures', 'PV_men.png'), width = 15, height = 10)

q6 <- panelview(death ~treatedpost , data=women, index=c("State","year"), pre.post= TRUE)

#ggsave(here('figures', 'PV_women.png'), width = 15, height = 10)

ggarrange(q5+ rremove("ylab") + rremove("xlab")+ggtitle("Men"), 
          q6+ rremove("ylab") + rremove("xlab")+ggtitle("Women"),
          ncol=1, nrow=2,common.legend = TRUE,
          cex.axis.x=20,
          cex.axis.y=20,
          cex.lab=20)
## $`1`

## 
## $`2`

## 
## $`3`

## 
## attr(,"class")
## [1] "list"      "ggarrange"
# ggsave(here('figures', 'efigure2_data_structure_men_women.png'), 
#        width = 40, height = 60, units = "cm")

5.10 eTable 3:Triple difference over time (interaction)

### Black or Hispanic vs White
ddd_black_hispanic_vs_white <- bind_rows(black_by_period %>% mutate(group="black"),
                hispanic_by_period %>% mutate(group="hispanic"),
                white_by_period %>% mutate(group="white")) %>% 
  dplyr::select(time, group, ATT, CI.lower, CI.upper, se=S.E.) %>% 
  pivot_wider(., names_from=group, 
              values_from=c(ATT, CI.lower, CI.upper, se)) %>% 
  summarise(mean_black=ATT_black-ATT_white,
            se_black=sqrt(se_black^2 + se_white^2),
            lower_black=round((mean_black-qnorm(p=1-(0.05)/2, mean=0, sd=1)*se_black), 2),
            upper_black=round((mean_black+qnorm(p=1-(0.05)/2, mean=0, sd=1)*se_black), 2),
            pvalue_black=round(pnorm(0,mean=abs(mean_black),sd=se_black) + 
              (1 - pnorm(0,mean=-abs(mean_black),sd=se_black)),2),
            mean_hispanic=ATT_hispanic-ATT_white,
            se_hispanic=sqrt(se_hispanic^2 + se_white^2),
            lower_hispanic=round((mean_hispanic-qnorm(p=1-(0.05)/2, mean=0, sd=1)*se_hispanic), 2),
            upper_hispanic=round((mean_hispanic+qnorm(p=1-(0.05)/2, mean=0, sd=1)*se_hispanic), 2),
            pvalue_hispanic=round(pnorm(0,mean=abs(mean_hispanic),sd=se_hispanic) + 
              (1 - pnorm(0,mean=-abs(mean_hispanic),sd=se_hispanic)),2),
            time=time) %>% 
  relocate(time) %>% pivot_longer(cols=-time,
               names_to = c(".value", "race"),
               names_sep="_") %>% 
  arrange(race)


### Women vs men
ddd_women_vs_men <- bind_rows(men_by_period %>% mutate(group="men"),
                women_by_period %>% mutate(group="women")) %>% 
  dplyr::select(time, group, ATT, CI.lower, CI.upper, se=S.E.) %>% 
  pivot_wider(., names_from=group, 
              values_from=c(ATT, CI.lower, CI.upper, se)) %>% 
  summarise(mean=ATT_women-ATT_men,
            se=sqrt(se_men^2 + se_women^2),
            lower=round((mean-qnorm(p=1-(0.05)/2, mean=0, sd=1)*se), 2),
            upper=round((mean+qnorm(p=1-(0.05)/2, mean=0, sd=1)*se), 2),
            pvalue=round(pnorm(0,mean=abs(mean),sd=se) + 
              (1 - pnorm(0,mean=-abs(mean),sd=se)),2),
            time=time) %>% 
  relocate(time)

ddd_table <- bind_rows(ddd_black_hispanic_vs_white, ddd_women_vs_men) %>% 
  rename(group = race) %>% 
  mutate(group = case_when(group=="black"~"Black vs White",
                           group=="hispanic"~"Hispanic vs White",
                           TRUE~"Women vs men")) %>% 
  mutate(`Adjusted Difference in Mean Difference (95%CI)`= simulr::ci_fct(
    pe=mean,
    low=lower,
    high=upper,
    digit=2,
    multiplier = 1,
    scale = "")) %>% 
  dplyr::select(Group= group, Time=time, `Adjusted Difference in Mean Difference (95%CI)`, `P-value`= pvalue)

# ddd_table %>% flextable::flextable() %>% 
#   flextable::save_as_docx(.,path=here("tables", "etable3.docx"))

ddd_table %>% 
  knitr::kable() %>%
  kableExtra::kable_styling(bootstrap_options = c("striped", "bordered"),
                full_width = F)
Group Time Adjusted Difference in Mean Difference (95%CI) P-value
Black vs White 0 9.97 (-4.92, 24.86) 0.19
Black vs White 1 9.88 (-6.97, 26.74) 0.25
Black vs White 2 2.61 (-21.46, 26.68) 0.83
Black vs White 3 7.05 (-20.84, 34.95) 0.62
Black vs White 4 8.33 (-16.89, 33.55) 0.52
Black vs White 5 11.89 (-12.25, 36.04) 0.33
Black vs White 6 16.16 (-7.19, 39.51) 0.17
Hispanic vs White 0 -4.01 (-16.60, 8.59) 0.53
Hispanic vs White 1 -4.28 (-22.71, 14.15) 0.65
Hispanic vs White 2 -2.68 (-25.43, 20.07) 0.82
Hispanic vs White 3 -4.03 (-31.20, 23.15) 0.77
Hispanic vs White 4 1.91 (-24.40, 28.22) 0.89
Hispanic vs White 5 -6.33 (-30.93, 18.27) 0.61
Hispanic vs White 6 -2.71 (-28.47, 23.05) 0.84
Women vs men 0 -2.90 (-10.76, 4.97) 0.47
Women vs men 1 -1.48 (-12.54, 9.57) 0.79
Women vs men 2 1.97 (-10.95, 14.89) 0.77
Women vs men 3 -6.86 (-20.18, 6.46) 0.31
Women vs men 4 1.71 (-13.08, 16.50) 0.82
Women vs men 5 -4.84 (-20.77, 11.10) 0.55
Women vs men 6 0.89 (-15.79, 17.57) 0.92